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Abstract 

A data compression algorithm involving vector quantization (VQ) and the discrete 
wavelet transform (DWT) is applied to two different types of multidimensional digital 
earth-science data. The algorithm (WVQ) is optimized for each particular application 
through an optimization procedure that assigns VQ parameters to the wavelet trans- 
form subbands subject to constraints on compression ratio and encoding complexity. 

Preliminary results of compressing global ocean model data generated on a Think- 
ing Machines CM-200 supercomputer are presented. The WVQ scheme is used in 
both a predictive and nonpredictive mode. Parameters generated by the optimiza- 
tion algorithm are reported, as are signal- to-noise ratio (SNR) measurements of actual 
quantized data. The problem of extrapolating hydrodynamic variables across the con- 
tinental landmasses in order to compute the DWT on a rectangular grid is discussed. 

Results are also presented for compressing Landsat TM 7-band data using the WVQ 
scheme.The formulation of the optimization problem is presented along with SNR 
measurements of actual quantized data. Postprocessing applications are considered 
in which the seven spectral bands are clustered into 256 clusters using a A’- means 
algorithm and analyzed using the Los Alamos multispectral data analysis program, 
SPECTRUM, both before and after being compressed using the WVQ program. 


I. Introduction. 

This work describes the application of an image compression algorithm involving the 
discrete wavelet transform and vector quantization to two problems involving earth science 
data. The coding of outputs of supercomputer-generated global climate model (GCM) ocean 
simulations and Landsat Thematic Mapper (TM) multispectral imagery is investigated. The 
compression algorithm has its origins in the coding of gray-scale imagery [1 ,2]. A set of vector 
quantizers is designed (one for each subband in the wavelet decomposition) with parameters 
selected from the solution of an optimization problem that is formulated to minimize quanti- 
zation distortion with constraints on the overall bit rate and encoding complexity. Although 
both data types considered in this work are of dimensionality higher than two, we restrict 
the discussion to coding implementations based on 2-D transforms. 
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Compression of the GCM data is approached by both a straightforward two-dimensional 
extension of the earlier algorithm and a predictive scheme in which two-dimensional pre- 
diction residuals are coded. The Landsat images are coded by modifying the bit allocation 
algorithm to allocate coder resources simultaneously among all of the spectral components. 
For all scenarios, measurements of quantization distortion are presented as functions of 
compression ratio and encoding complexity, thus revealing tradeoffs involved in the system 
design. 


II. Multidimensional Wavelet Transform-Vector Quantization. 

The data-coding technique used in this work, known as the wavelet-vector quantization 
(WVQ) algorithm, is based on vector quantization of the subbands resulting from a discrete 
wavelet transform (DWT) decomposition of the data signal. For signals in two or more 
dimensions, the transform used is based on product filter banks (i.e., tensor products of 
one-dimensional DWT filters). A d-dimensional signal transformed in this manner with a 
two-channel filter bank yields 2 d subbands, any of which can be cascaded back through 
the filter bank to produce a multirate decomposition of the original signal. Although we are 
currently working on three-dimensional wavelet transforms for use with the three-dimensional 
climate model data under investigation, the DWT results presented here are restricted to 
the case of two-dimensional data fields. 

Single-level 2-D DWT analysis and synthesis filter banks are depicted in Figures 1 and 
2. The analysis filters (Hi) and synthesis filters (F*) used in this paper are biorthogonal 
linear phase FIR wavelet filters constructed in [3, 4]. Note the use of binary subscripts on 
the subbands, a*j, to indicate the filters applied to the rows and columns of the signal, x. 
Signals obtained from sampling smooth, continuous data fields usually have most of their 
energy (or variance) concentrated in the low-frequency part of the spectrum, so it is usually 
most efficient to cascade only the lowpass-lowpass filtered subband, a 0 cu back through the 
analysis bank in Figure 1; the resulting subband is denoted a 0 o,oo> or a )0 o f° r short. This 
cascade is typically carried down four levels (or so), to the a m oo band, which will then contain 
a large portion of the signal energy concentrated in a heavily downsampled signal component. 
Note that the downsample factor for an Z^-level subband in a d-dimensional scheme using 
an M-channel coder is ra t = M ld , so, e.g., subband a nt00 in the 2-D decomposition has been 
downsampled by m t = 256-to-l. 

A further consideration when applying a DWT to finite-duration signals, like the rows or 
columns in a digital image, is the handling of boundary conditions. The most straightforward 
way of dealing with signal boundaries is to regard the signal as a single period of an infinite, 
periodic input and apply the DWT filter bank by circular convolution and downsampling. 
This has the disadvantages, however, of introducing a spurious jump discontinuity when the 
data isn’t inherently periodic and of constraining the signal length (i.e., its “period”) to 
be divisible by the downsample factor. For a four-level decomposition using a two-channel 
filter bank, for instance, this means the input length, iY 0 , must be divisible by 16. Both 
of these problems can be avoided by using symmetric extrapolation techniques to extend 
finite-duration inputs; moreover, this can be done with no increase in the memory allocation 
needed to transform or store the input signal [5, 6]. 

The design of vector quantizers for the subbands in a DWT decomposition is based on 
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the assumption of exponential VQ rate-distortion characteristics, 

D i (k i ,r i ) = / 3 i (k t )e-^ r ' (1) 

D^ki, 77) is the distortion (mean-square error) between the original and quantized data in 
the i th subband for a bit rate of r t bits per pixel (bpp) and a vector dimension of Ar,. /?,(&,) 
and 7 ,(&,) are constants that depend on k, and the probability density function of the data 
vectors. The motivation for this assumption is based on theoretical VQ rate-distortion 
modelling [ 7 ] and confirmed by empirical data; values for the constants and 7 t (k t ) are 

determined from a set of training data. 

In the case of an orthogonal subband decomposition, the overall distortion can be ex- 
pressed as a weighted sum of the distortion in each subband, 


£> = £ — D t (k t ,r t ) , 

T m - 

where rn t is the downsample factor, the ratio of the number of samples in the original to 
the number in the i th subband. Since the DWT conserves the number of data samples, m, 
satisfies the identity = 1 - By (1), the overall distortion is 




Formula (2) is customarily used as a distortion measure with nonorthogonal transforms, too, 
although it no longer coincides exactly with overall mean-square error. The bit-allocation 
problem for quantizer design involves using nonlinear optimization techniques to compute 
the bit rates, 77, and dimensions, k t , that minimize (2), subject to constraints on overall bit 
rate and encoder complexity. 

For a target overall bit rate of R bpp, the constraint on subband bit rates is 



< R • 



If subband vector dimensions, Q, are to be optimized, an additional constraint besides 


( 3 ) is necessary to obtain a well-posed optimization problem for VQ bit allocation. The 
encoder complexity constraint used here is an upper bound, Q , on the computational cost of 
performing exhaustive nearest-neighbor searches of k , -dimensional VQ codebooks containing 
iV, — 2 k,r ' codevectors: 


Y — 2 k ' T 'a < Q 

tlli 



The parameter a is a constant corresponding to the arithmetic cost of performing two ad- 
ditions and one multiplication per pixel. With the additional constraints 77 > 0 we obtain a 
convex nonlinear optimization problem to solve for the 77; the k, are optimized by a heuris- 
tic search procedure. Once optimal bit rates and vector dimensions are computed, optimal 
VQ codebooks are constructed from training data using the Linde- Buzo- Gray method [8, 9]. 
More details about the WVQ algorithm are given in [ 10 , 2 , 1, 11], 
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III. Application to Ocean Model Data. 

This section describes the use of the WVQ algorithm on synthetic data generated by 
a Bryan- Cox- Semtner global ocean circulation model running on the Connection Machine 
CM-200 at the Los Alamos Advanced Computing Laboratory (ACL) [12, 13]. The model 
is computed on a 320 x 768 grid at 20 depth levels; boundary conditions are given on a 
three-dimensional bottom topography with 80 islands. The data used in the compression 
experiments was the surface temperature field (no depth components), taken at three-day 
intervals over a decade's worth of simulation. Time-frames from the first year of the simula- 
tion were used for training data, and the resulting WVQ algorithm was then tested on frames 
from the last year of the simulation, i.e., on data similar but not identical to the training 
data. We feel this is a valid test since it is similar to the manner in which the algorithm will 
be used in practice. 

The two-dimensional data frames were transformed with a four-level octave-scaled DWT 
decomposition. Since the model is periodic in the east-west direction, periodic boundary 
conditions were used along parallels of latitude (“rows”). However, due to the lack of con- 
tinuity between the north and south edges of the grid, symmetric (i.e., reflected) boundary 
conditions were used along meridians of longitude (“columns”). Because the DW T is most 
easily computed on a rectangular grid, the temperature data was extended across the con- 
tinental landmasses before transforming. A simple approach like zero-padding of the data 
would be undesirable because it would induce a large jump discontinuity around the coast- 
lines; this would show up as added variance in the highpass-filtered DW 1 subbands and 
would therefore reduce the compressibility of the high-frequency signal components in the 
transform domain. For this reason we used a continuous extension of the data given by linear 
extrapolation from coast to coast along parallels of latitude. This still leaves a corner at 
the coastlines in the extended data; since the initial data field is extremely smooth, this 
corner results in a slight increase in energy in highpass-filtered subbands. It is not yet clear 
whether this added variance is significant alongside the variance naturally present in the 
data. We are currently looking into using smoother two-dimensional extrapolation schemes 
for this task. 

Two different approaches were taken to quantizing the time-series data generated in the 
simulation: nonpredictive and predictive coding. In the nonpredictive scheme, each frame 
is treated as a separate image and compressed accordingly using the W 7 VQ method. In 
predictive coding, a prediction of each frame is made based on past frames and subtracted 
from the current frame, resulting in a two-dimensional residual image, which is compressed 
and stored. The image sequence is decoded from the first frame in the sequence and the 
residuals. We used a simple first-order predictor in this scheme; i.e., the prediction of a given 
frame is just equal to the quantized value of the previous frame. Block diagrams for the 
transmitter and receiver in this predictive encoding/decoding system are given in Figures 3 
and 4. The experiment assumed that the first frame in the sequence is transmitted with 
nonpredictive quantization, and the compression ratios reported are those of the subsequent 
residuals. 

For both the nonpredictive and predictive schemes, WVQ coders were designed for bit 
rates, R, ranging from 2.0 to 0.25 bpp. Since the original data was 32 bpp, the corresponding 
compression ratios range from 16:1 to 128:1. Encoding complexities, Q , were varied between 
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Figure 3: Predictive Transmitter. 




Predictor 


Figure 4: Predictive Receiver. 


16a and 64a. The optimal codebook sizes and vector dimensions for each combination of R 
and Q were computed by the WVQ design algorithm described in Section II. Bit allocation 
results for the 13 subbands are presented in Tables I and II for R = 0.5 and 0.25 bpp in 
terms of vector dimensions, and codebook sizes, N t . Note that as the bit rate decreases, 
the highest frequency subbands are quantized more heavily or discarded altogether and 
remaining high-frequency subband vector dimensions typically increase. Vector dimensions 
also increase as the upper bound on complexity increases. Bit allocation for the residual 
subbands in Table II is similar to that for the nonpredictive scheme, although much less of 
the quantizer resources (he., far fewer bits) are allocated to the lower frequency subbands in 
the predictive scheme. This means that the first-order predictor effectively predicts the low- 
wavenumber modes of the model, indicating that these modes are evolving slowly compared 
to the sampling rate for archiving data. 

Quantizer performance is measured in terms of signal- to-noise ratio (SNR), 

2 

SNR = 101og 10 (dB) , 

a i 

where a* is the signal variance and a\ is the quantization error variance. The average SNR 
in dB for the test data is shown in Figures 5 and 6 for various values of R and Q. This 
diagram illustrates the various tradeoffs involved in the selection of R and Q. At a given bit 
rate, /?, note that a higher SNR is possible using an encoder with higher complexity, Q ; i.e., 
higher subband vector dimensions, k t , and correspondingly larger codebook sizes, TV,. Note 
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Table I: Vector Dimension and Codebook Size Assignments ( k,N ) for R = 0.5 bpp and 
R — 0.25 bpp, Nonpredictive Coding. 



Q = 64 a 

R = 0.5 
Q = 32a 

Q = 16a 

Q = 64a 

R = 0.25 
Q = 32a 

Q = 16a 

Subband ,,,00 

(1,511) 

(1,493) 

(1,308) 

(1,251) 

(1,251) 

(1,251) 

Subband ,,,01 

(1,132) 

(1,119) 

(1,91) 

(1,32) 

(1,31) 

(1,32) 

Subband ,,,10 

(1,93) 

(1,103) 

(1,75) 

(1,15) 

(1,17) 

(1,18) 

Subband ,,,11 

(1,48) 

(1,58) 

(1,41) 

(1,8) 

(1,7) 

(1,7) 

Subband ,,01 

(4,723) 

(2,247) 

(2,133) 

(4,377) 

(4,382) 

(4,295) 

Subband ,,10 

(4,261) 

(2,61) 

(2,41) 

(4,16) 

(4,16) 

(4,17) 

Subband ,,11 

(4,134) 

(4,76) 

(2,25) 

(4,4) 

(4,4) 

(4,4) 

Subband ,01 

(8,373) 

(8,177) 

(4,68) 

(8,176) 

(8,176) 

(8,132) 

Subband ,10 

(8,19) 

(8,13) 

(8,8) 

— 

— 

— 

Subband ,1 1 

— 

— 

(8,2) 

— 

— 

— 

Subband 01 

(8,74) 

(8,42) 

(8,22) 

(8,4) 

(8,4) 

(8,4) 

Subband 10 

— 

— 

— 

— 

— 

— 

Subband 1 1 

— 

— 

— 

— 

— 

— 


Table II: Vector Dimension and Codebook Size Assignments ( k,N ) for R = 0.5 bpp and 
R = 0.25 bpp, Predictive Coding. 



Q = 64 a 

R = 0.5 
Q = 32a 

Q ~ 16a 

Q = 64a 

R = 0.25 
Q — 32 a 

Q = 16a 

Subband ,,,00 

(1,14) 

(1,10) 

(1,16) 

(1,2) 

(1,2) 

(1,4) 

Subband ,,,01 

(1,21) 

(1,15) 

(1,23) 

( 1 , 3 ) 

(1,4) 

(1,6) 

Subband ,,,10 

(1,26) 

(1,20) 

(1,28) 

(1,5) 

(1,6) 

( 1 , 9 ) 

Subband ,„11 

(1,35) 

(1,28) 

(1,35) 

(1,8) 

(1,10) 

( 1 , 14 ) 

Subband ,,01 

(4,383) 

(4,235) 

(2,84) 

(4,207) 

(4,198) 

(4,120) 

Subband ,,10 

(4,423) 

(2,125) 

(2,86) 

(4,293) 

(4,257) 

(4,140) 

Subband ,,11 

(4,305) 

(4,185) 

(2,64) 

(4,93) 

(4,109) 

(4,87) 

Subband ,01 

(8,276) 

(4,126) 

(4,62) 

(8,461) 

(8,209) 

(8,78) 

Subband ,10 

(8,239) 

(8,138) 

(8,54) 

(8,206) 

(8,145) 

(8,65) 

Subband ,11 

(8,97) 

(8,51) 

(8,24) 

— 

(8,3) 

(8,12) 

Subband 01 

(8,30) 

(8,12) 

(8,11) 

— 

— 

— 

Subband 10 

— 

— 

— 

— 

— 

— 

Subband 11 

— 

— 

— 

— 

— 

— 
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Figure 5: SNR Measurements of Quantized Temperature Data, Nonpredictive. 
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Figure 6: SNR Measurements of Quantized Temperature Data, Predictive Coding 




that the complexity can he increased in this manner while the subband bit rates, r*, remain 
unchanged since 

log 2 N i 


Ti ~ 


h 


Increasing the encoding complexity results in an encoder with a more time-consuming code- 
book lookup but does not affect decoder performance. As the bit rate increases, we see from 
Figure 5 that the gain in SNR achieved by increasing the encoding complexity becomes more 
significant. For a fixed R and Q , comparison of Figures 5 and 6 shows that the predictive 
scheme achieves a gain on the order of 1-4 dB over nonpredictive coding. The improvement 
in coding gain is more pronounced at lower bit rates and higher complexities, since the resid- 
ual coding scheme is better able to exploit higher limits on encoding complexity at low bit 
rates than the nonpredictive scheme. 


IV. Application to Multispectral Data. 


This section discusses the application of WVQ to the compression of multispectral im- 
agery. Since each spectral band is a separate monochromatic image, the approach is to code 
each of the bands by two-dimensional WVQ using symmetric boundary conditions. The 
bit-allocation is performed for the various spectral components simultaneously and hence 
the coding of each spectral component is not viewed as a separate two-dimensional problem. 

The multispectral problem requires a modification to the WVQ design algorithm dis- 
cussed in Section II since the rate and complexity are expressed in terms of multidimensional 
pixels. For the case of L spectral components the system design procedure entails minimizing 


D = W—m)e-^ )T ' 

Lf ; 771 i 


(5) 


over the k{ and i\ subject to 


t£-<* 

L i mi 

IE— V k 'a<Q 

L 777 ; 


T { > 0 
k» 6; A; 


( 6 ) 

(7) 

( 8 ) 
(9) 


where K t denotes a prespecified set from which k i must be selected. The optimization 
is performed over all of the two-dimensional subbands generated from all of the spectral 
components. 

Multispectral image WVQ was considered for the application of compressing Landsat 
Thematic Mapper (TM) data. Such data consist of seven 8-bit spectral bands (three visible, 
three infrared, and one thermal) at a ground sample distance of 28.5 meters. Four data sets 
were used in training the coder: Albuquerque, NM (2984 x 3356); Cairo, Egypt (2945 x 3320); 
Los Alamos, NM (2984 x 3254); and Mexico City, Mexico (5965 x 6967). The performance 
of the coder was evaluated in terms of results obtained by compressing a (2976 x 3552) scene 
from the Moscow, Russia, area containing both urban and agricultural areas. The resulting 
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Table III: RMSE Quantizer Performance as a Function of Compression Ratio and Complex- 
ity. 



Q = 8a 

Q = 16a 

Q = 32a 

Q = 64a 

16:1 

3.04 

2.35 

1.92 

1.70 

32:1 

3.06 

2.63 

2.43 

2.26 

64:1 

3.85 

3.51 

3.32 

3.24 

128:1 

4.76 

4.71 

4.58 

4.56 


root mean-square error (RMSE) for sixteen combinations of bit rate and encoding complexity 
are tabulated in Table III. The compression ratios reported are relative to 56 bits/pixel in 
the original data and assume that the bit rates satisfy 


t£- = * 

1 mi 


with L — 1 


( 10 ) 


The additional gain available from entropy and run-length coding is not included. 

It is interesting to compare these results to those obtained by another wavelet-based 
compression technique. In [14] Landsat TM images were compressed via a subband de- 
composition of each spectral component by a 7-tap nonperfect reconstruction filter bank; 
each subband was coded with uniform scalar quantization followed by Huffman and zero- 
run-length coding. The experiment was repeated with an image-dependent Karhunen-Loeve 
transform (KLT) in the interband direction, which provided noticeable coding gain at the 
expense of computational complexity. The WVQ RMSE results depicted in Table III appear 
to lie between these two previous approaches, although any comparisons must be qualified by 
the fact that the numerical results in these two papers were obtained from different imagery 
(Kuwaiti oil fields in the case of [14]). For instance, Table III shows that 32:1 compression 
(1.75 bpp) with a complexity of Q = 64a yields an average RMSE per band of 2.26, or a 
little over 2 bits of error. The closest comparable value for non-KLT coding in [14] is a MSE 
of 40.02 at 2.51 bpp; dividing the MSE by 7 and taking a square root gives an average RMSE 
per band of 2.39, which is a slightly greater error at a higher bit rate than our result. With 
interband KLT coding, [14] reports a MSE of 25.11, or an average RMSE of 1.89, at 1.55 
bpp; this is a lower distortion at a lower bit rate than our result. We are currently currently 
working on incorporating interband KLT coding with the WVQ compression method. 

The motivation for our investigation of TM data compression is the need to store and 
process large amounts of data for postprocessing applications. Using the software package 
SPECTRUM [15], developed by Los Alamos National Laboratory and the University of New 
Mexico, we are able to use a desktop workstation running Unix and X-windows to analyze 
and categorize multispectral data that has been clustered into 256 clusters using a variant of 
the fc-means algorithm. SPECTRUM can manipulate the color map for the computer display 
using any transformation of the clustered data, and can display cluster position as a two- 
dimensional scatter plot. Using these features, users are able to categorize data by selecting 
areas with a known type of land cover, causing all associated pixels in the image to be given 
the same pseudocolor representation. Of great interest to us is the robustness of SPECTRUM 
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data clustering when applied to data that has first been compressed by the WVQ algorithm. 
While visual quality of pseudocolor visualizations remains good after compression by as 
much as 32:1, it remains to be determined how much quantization distortion SPECTRUM 
can tolerate for tasks like Level 1 Land Use Categorization. We are attempting to establish 
quantitative distortion criteria based on the analysis of classification error presented in [15], 
which is based on computing levels of confidence for classifications done on clustered data. 
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